#include <stdio.h>
#include <assert.h>

#define cudaCheckErrors(msg) __getLastCudaError(msg, __FILE__, __LINE__)

inline void __getLastCudaError(const char *errorMessage, const char *file,
                               const int line) {
  cudaError_t err = cudaGetLastError();

  if (cudaSuccess != err) {
    fprintf(stderr,
            "%s(%i) : CUDA error :"
            " %s : (%d) %s.\n",
            file, line, errorMessage, static_cast<int>(err),
            cudaGetErrorString(err));
    exit(EXIT_FAILURE);
  }
}



// Convenience function for checking CUDA runtime API results
// can be wrapped around any runtime API call. No-op in release builds.
inline
cudaError_t checkCuda(cudaError_t result)
{
#if defined(DEBUG) || defined(_DEBUG)
  if (result != cudaSuccess) {
    fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
    assert(result == cudaSuccess);
  }
#endif
  return result;
}

const int BLOCK_SIZE = 256;
const int NUM_REPS   = 4;

// Check errors and print GB/s
void postprocess(long int n, float ms,int stride){
  printf("                  %4d %20.2f\n", stride,2 * n * sizeof(float) * 1e-6 * NUM_REPS / ms );
}


// simple copy kernel
__global__ void copy(float *odata, const float *idata,int stride){
  int i = blockIdx.x * blockDim.x + stride * threadIdx.x;
  odata[i] = idata[i];
}


int main(int argc, char **argv){
  const long int nx     = 1024L * 1024 * 256;
  const size_t mem_size = nx*sizeof(float);

  dim3 dimGrid(nx/BLOCK_SIZE,1,1);
  dim3 dimBlock(BLOCK_SIZE,1,1);
  int stride = 1;

  int devId = 0;
  if (argc > 1) devId = atoi(argv[1]);

  cudaDeviceProp prop;
  checkCuda( cudaGetDeviceProperties(&prop, devId));
  printf("\nDevice : %s\n", prop.name);
  printf("Array size: %d, Block size: %d\n",nx,BLOCK_SIZE);
  printf("dimGrid: %d %d %d. dimBlock: %d %d %d\n",
         dimGrid.x, dimGrid.y, dimGrid.z, dimBlock.x, dimBlock.y, dimBlock.z);
  
  checkCuda( cudaSetDevice(devId) );

  float *h_idata = (float*)malloc(mem_size);
  float *h_cdata = (float*)malloc(mem_size);
  float *d_idata, *d_cdata;
  checkCuda( cudaMalloc(&d_idata, mem_size) );
  checkCuda( cudaMalloc(&d_cdata, mem_size) );

  // check parameters and calculate execution configuration
  if (nx % BLOCK_SIZE) {
    printf("nx must be a multiple of BLOCK_SIZE\n");
    goto error_exit;
  }

  // host
  for (int i = 0; i < nx; i++) h_idata[i] = sin(1.0e-8*i);
  
  // device
  checkCuda( cudaMemcpy(d_idata, h_idata, mem_size, cudaMemcpyHostToDevice) );
  
  // events for timing
  cudaEvent_t startEvent, stopEvent;
  checkCuda( cudaEventCreate(&startEvent) );
  checkCuda( cudaEventCreate(&stopEvent) );
  float ms;

  // ------------
  // time kernels
  // ------------
  printf("%25s%25s\n", "Stride", "Bandwidth (GB/s)");
  
  // ----
  // copy 
  // ----
  //printf("%25s", "copy");
  checkCuda( cudaMemset(d_cdata, 0, mem_size) );
  // warm up
  copy<<<dimGrid,dimBlock>>>(d_cdata,d_idata,stride);
  cudaCheckErrors("Warm up");
  for (stride=1;stride<256;stride++){
    checkCuda( cudaEventRecord(startEvent, 0) );
    checkCuda( cudaEventSynchronize(startEvent) );
    for (int i = 0; i < NUM_REPS; i++)
      copy<<<nx/BLOCK_SIZE-stride,BLOCK_SIZE>>>(d_cdata, d_idata,stride);
    cudaCheckErrors("Copies");
    checkCuda( cudaEventRecord(stopEvent, 0) );
    checkCuda( cudaEventSynchronize(stopEvent) );
    checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
    checkCuda( cudaMemcpy(h_cdata, d_cdata, mem_size, cudaMemcpyDeviceToHost) );
    postprocess(nx-stride*BLOCK_SIZE,ms,stride);
  }
 


error_exit:
  // cleanup
  checkCuda( cudaEventDestroy(startEvent) );
  checkCuda( cudaEventDestroy(stopEvent) );
  checkCuda( cudaFree(d_cdata) );
  checkCuda( cudaFree(d_idata) );
  free(h_idata);
  free(h_cdata);
}
